{
 "cells": [
  {
   "cell_type": "code",
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:55.561411Z",
     "start_time": "2025-08-09T02:52:55.557374Z"
    }
   },
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import arange\n",
    "\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题"
   ],
   "execution_count": 66,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:55.907913Z",
     "start_time": "2025-08-09T02:52:55.904277Z"
    }
   },
   "cell_type": "code",
   "source": [
    "param = np.array([\n",
    "    [30,50,6,18.193],\n",
    "    [30,70,8,9.12],\n",
    "    [30,90,10,9.276],\n",
    "    [40,50,6,30.243],\n",
    "    [40,70,8,23.856],\n",
    "    [40,90,10,35.673],\n",
    "    [50,50,6,6.586],\n",
    "    [50,70,8,6.586],\n",
    "    [50,90,10,8.323],\n",
    "])\n"
   ],
   "id": "df285f9baf5dba76",
   "execution_count": 67,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:56.171733Z",
     "start_time": "2025-08-09T02:52:56.168937Z"
    }
   },
   "cell_type": "code",
   "source": "",
   "id": "2b850109c40f3fb3",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:56.353213Z",
     "start_time": "2025-08-09T02:52:56.343650Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def f(T, H, SC, klist, khlist):\n",
    "    k_1 = klist[0]\n",
    "    k_3 = klist[1]\n",
    "    k_4 = klist[2]\n",
    "    k_5 = klist[3]\n",
    "    k_6 = klist[4]\n",
    "\n",
    "    S_C = khlist[0]\n",
    "    S_S = khlist[1]\n",
    "    k_H = khlist[2]\n",
    "    # 变量\n",
    "    solid_content = SC / 100  # 固含量\n",
    "    T = T + 273.15  # 温度，单位为开尔文\n",
    "\n",
    "    # 常量\n",
    "    # 常量\n",
    "    C_0 = 24e-3  # 任意的比例系数，用于确定DMF质量与时间的关系\n",
    "    m_D_0 = 24e-3  # 初始时刻DMF的质量，单位为kg\n",
    "    m_S = 6e-3  # 环丁砜总的质量，单位为g。无论是否溶解，环丁砜不会从体系中消失，总的质量是不变的\n",
    "    m_C = (m_D_0 + m_S) * solid_content  # 初始时刻醋酸纤维素的质量，单位为kg，与初始固含量有关\n",
    "    ro_D = 0.948e3  # DMF的密度，单位为kg/m^3\n",
    "    ro_S = 1.261e3  # 环丁砜的密度，单位为kg/m^3\n",
    "    ro_C = 1.3e3  # 醋酸纤维素的密度，单位为kg/m^3\n",
    "    V = m_D_0 / ro_D + m_S / ro_S + m_C / ro_C  # 总体积，单位为m^3，认为不变\n",
    "    A = 6.09451\n",
    "    B = 2725.96\n",
    "    C = 28.209 # 由温度计算饱和汽压时的三个常数\n",
    "    A_0 = -k_1 * (10 ** (A - B / (T + C)))*133.322 * (1 - k_H * H / 100) / ro_D / V\n",
    "    n_correction_factor = 6  # 斯托克斯 - 爱因斯坦方程中的修正系数，对于大分子乃至宏观颗粒是6，对于小分子则体积越小，该值也会越小\n",
    "    viscosity_rate = 0.92e-3  # 在20℃下DMF粘度与水粘度的比值\n",
    "    k = 1.380649e-23  # 玻尔兹曼常数，单位为J/K\n",
    "    r = 0.3413e-9  # 环丁砜的分子半径\n",
    "    N_A = 6.02214076e23  # 阿伏伽德罗常量，单位为mol^-1\n",
    "\n",
    "    # 函数\n",
    "    # DMF在蒸发的过程中，其质量随时间变化的函数\n",
    "    def m_D(t):\n",
    "        return C_0 * np.exp(A_0 * t)\n",
    "\n",
    "\n",
    "    # 未能溶解而析出成液滴的环丁砜的质量随时间变化的函数\n",
    "    def m_S_out(t):\n",
    "        return max(0, m_S - m_D(t) * S_S)\n",
    "\n",
    "\n",
    "    # 未能溶解而析出成固体的醋酸纤维素的质量随时间变化的函数\n",
    "    def m_C_out(t):\n",
    "        return max(0, m_C - m_D(t) * S_C)\n",
    "\n",
    "\n",
    "    # 析出的醋酸纤维素的体积占比随时间变化的函数\n",
    "    def phi_C_out(t):\n",
    "        return m_C_out(t) / ro_C / V\n",
    "\n",
    "\n",
    "    # 由水在30、40、50（℃）下的粘度和DMF与水粘度的关系计算不同温度下DMF的粘度，单位为mPa*s\n",
    "    def eta_0(temp):\n",
    "        if temp == 30+273.15:\n",
    "            return 0.8007 * viscosity_rate\n",
    "        elif temp == 40+273.15:\n",
    "            return 0.6560 * viscosity_rate\n",
    "        elif temp == 50+273.15:\n",
    "            return 0.5494 * viscosity_rate\n",
    "\n",
    "\n",
    "    # 不同温度下混合溶液总的粘度随时间变化的函数\n",
    "    def eta(temp, t):\n",
    "        return eta_0(temp) * (1 + 2.5 * phi_C_out(t))\n",
    "\n",
    "\n",
    "    # 扩散系数。假定一个小液滴内仅有一个分子，且分子为球形\n",
    "    def D(temp, t):\n",
    "        return k * temp / (n_correction_factor * np.pi * r * eta(temp, t))\n",
    "\n",
    "\n",
    "    # 液滴运动的平均速度\n",
    "    def v(temp, t):\n",
    "        return k_3 * (D(temp, t) ** 0.5)\n",
    "\n",
    "\n",
    "    # 析出的环丁砜的分子数密度随时间变化的函数\n",
    "    def n_molecular_number_density(t):\n",
    "        return m_S_out(t) / (ro_S*4/3*np.pi*r**3) / V\n",
    "\n",
    "\n",
    "    # 小液滴间的碰撞频率\n",
    "    def Z(temp, t):\n",
    "        return (2 ** 0.5) * n_molecular_number_density(t) * np.pi * (r ** 2) * v(temp, t)\n",
    "\n",
    "    tf = 0\n",
    "    for t in np.arange(0, 500, 0.01):\n",
    "\n",
    "        if m_D(t) < 0.024*0.1:\n",
    "            tf = t\n",
    "            break\n",
    "    # print(\"tf：\",tf)\n",
    "\n",
    "    dy = lambda m_s, t: k_4 * Z(T, t) + k_5 * v(T, t) * (m_S_out(t) - m_s) / V * m_s\n",
    "    t = arange(1, tf, 0.01)\n",
    "    sol = odeint(dy, 0, t)\n",
    "\n",
    "    if len(sol.T[0]) == 0:\n",
    "        re = 10000\n",
    "    else:\n",
    "        re = k_6*sol.T[0][-1]\n",
    "    # print(sol.T[0][-1])\n",
    "    return re\n",
    "\n",
    "print(f(50,90,10,[1.79848413, 3.49918172, 0.55301977, 7.33290613,1],[0.9,0.9,0.2]))\n"
   ],
   "id": "4598c346b7688cc5",
   "execution_count": 68,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:56.550037Z",
     "start_time": "2025-08-09T02:52:56.545105Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def pred(klist,khlist):\n",
    "    predictions = []\n",
    "    for i in range(param.shape[0]):\n",
    "        T = param[i][0]\n",
    "        H = param[i][1]\n",
    "        SC = param[i][2]\n",
    "        re = f(T,H,SC,klist,khlist)\n",
    "        predictions.append(re)\n",
    "    return predictions\n",
    "print(pred([1,1,1,1,1],[0.9,0.9,0.2]))"
   ],
   "id": "1cf1e6dd5880b7b5",
   "execution_count": 69,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:56.737024Z",
     "start_time": "2025-08-09T02:52:56.733343Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# def target(T,H,SC):\n",
    "#     tar = 0\n",
    "#     for i in range(9):\n",
    "#         if param[i,0]==T and param[i,1]==H and param[i,2]==SC:\n",
    "#             tar = param[i,3]\n",
    "#     return tar\n",
    "# print(target(30,50,8))\n",
    "def target():\n",
    "    targets = []\n",
    "    for i in range(param.shape[0]):\n",
    "        t = param[i][-1]\n",
    "        targets.append(t)\n",
    "    return targets\n",
    "print(target())"
   ],
   "id": "7fec16d32cdfff96",
   "execution_count": 70,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:57.141619Z",
     "start_time": "2025-08-09T02:52:57.119412Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# def loss(pred, target):\n",
    "#     return np.mean((pred - target)**2)\n",
    "\n",
    "def loss(pred, target):\n",
    "    loss = 0\n",
    "    for i in range(param.shape[0]):\n",
    "        loss += np.linalg.norm(pred[i] - target[i])\n",
    "    return loss\n",
    "print(loss(pred([0.01646966, 0.03255624, 0.06004486, 0.07621205, 0.97623144],[0.85945502, 0.89168254, 0.39898241]),target()))"
   ],
   "id": "cba65906f0368f5a",
   "execution_count": 71,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:57.653382Z",
     "start_time": "2025-08-09T02:52:57.650024Z"
    }
   },
   "cell_type": "code",
   "source": [
    "kh_bound = np.array([\n",
    "    [0.7,1],\n",
    "    [0.7,1],\n",
    "    [0,0.8],\n",
    "])\n",
    "\n",
    "k_bound = np.array([\n",
    "    [0,1],\n",
    "    [0.001,1],\n",
    "    [0.001,1],\n",
    "    [0.001,1],\n",
    "    [0.001,1],\n",
    "])"
   ],
   "id": "ec0ff9c50f53d774",
   "execution_count": 72,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:58.076566Z",
     "start_time": "2025-08-09T02:52:58.072710Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def printre(klist,khlist):\n",
    "    predictions = []\n",
    "    t = target()\n",
    "    for i in range(param.shape[0]):\n",
    "        T = param[i][0]\n",
    "        H = param[i][1]\n",
    "        SC = param[i][2]\n",
    "        re = f(T,H,SC,klist,khlist)\n",
    "        tt = t[i]\n",
    "        print(f\"预测值{re},真实值{tt}\")"
   ],
   "id": "73dd088928c78e1c",
   "execution_count": 73,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:58.531644Z",
     "start_time": "2025-08-09T02:52:58.526110Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plot_SA_history():\n",
    "    plt.figure(figsize=(14, 10))\n",
    "\n",
    "    # 子图1：损失函数变化\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.semilogy(history['best_loss'], 'r-', label='Best Loss')\n",
    "    plt.semilogy(history['current_loss'], 'b--', alpha=0.5, label='Current Loss')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Loss')\n",
    "    plt.title('Loss变化曲线')\n",
    "    plt.legend()\n",
    "    plt.grid(True, which=\"both\", ls=\"--\")\n",
    "\n",
    "    # 子图2：温度衰减曲线\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.plot(history['temperature'], 'g-')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Temperature')\n",
    "    plt.title('温度下降曲线')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k参数演化\n",
    "    plt.subplot(2, 2, 3)\n",
    "    k_array = np.array(history['k_params'])\n",
    "    for i in range(5):\n",
    "        plt.plot(k_array[:, i], label=f'k_{i+1}')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k参数值')\n",
    "    plt.title('k参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图4：kh参数演化\n",
    "    plt.subplot(2, 2, 4)\n",
    "    kh_array = np.array(history['kh_params'])\n",
    "    for i in range(3):\n",
    "        plt.plot(kh_array[:, i], label=f'kh_{i+1}')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('物理参数值')\n",
    "    plt.title('物理参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "b268117688e30fd4",
   "execution_count": 74,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:52:59.347282Z",
     "start_time": "2025-08-09T02:52:59.343542Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 在SA函数外部初始化记录容器\n",
    "history = {\n",
    "    'temperature': [],\n",
    "    'best_loss': [],\n",
    "    'current_loss': [],\n",
    "    'k_params': [],\n",
    "    'kh_params': []\n",
    "}"
   ],
   "id": "73818db199139138",
   "execution_count": 75,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:56:43.085318Z",
     "start_time": "2025-08-09T02:56:30.628712Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def SA(t0,tf,alpha,iter):\n",
    "    global history\n",
    "    flag = 0\n",
    "    t = t0\n",
    "    k0 = np.array([0.01,0.1,0.02,0.1,0.1])\n",
    "    kh0 = np.array([0.9,0.9,0.3])\n",
    "    kc = k0\n",
    "    khc = kh0\n",
    "    lc = loss(pred(kc,khc),target())\n",
    "\n",
    "    kb = k0\n",
    "    lb = 1000\n",
    "    for i in range(iter):\n",
    "        kn = kc + np.random.normal(0,0.02,size=5)\n",
    "        kn = np.clip(kn,k_bound[:,0],k_bound[:,1])\n",
    "        khn = khc + np.random.uniform(-0.05,0.05,size=3)\n",
    "        khn = np.clip(khn,kh_bound[:,0],kh_bound[:,1])\n",
    "        pre = pred(kn,khn)\n",
    "        ln = loss(pre,target())\n",
    "\n",
    "        # for p in pre:\n",
    "        #     if p < 0.001:\n",
    "        #         flag = 1\n",
    "        # if flag==1:\n",
    "        #     continue\n",
    "\n",
    "        if ln < lc or np.random.rand() < np.exp(-(ln-lc)/t):\n",
    "            kc = kn\n",
    "            khc = khn\n",
    "            lc = ln\n",
    "            if lc < lb:\n",
    "                kb = kc\n",
    "                khb = khn\n",
    "                lb = lc\n",
    "        flag = 0\n",
    "        t *= alpha\n",
    "        if t < tf:\n",
    "            break\n",
    "        print(f\"iter{i}, lb{lb},kc{kc},khc{khc},lc{lc},\")\n",
    "        history['temperature'].append(t)\n",
    "        history['best_loss'].append(lb)\n",
    "        history['current_loss'].append(lc)\n",
    "        history['k_params'].append(kc.copy())\n",
    "        history['kh_params'].append(khc.copy())\n",
    "    return kb,khb\n",
    "re = SA(300,0.001,0.95,200)\n",
    "print(re)\n",
    "plot_SA_history()"
   ],
   "id": "36d7d95fa58cf227",
   "execution_count": 79,
   "outputs": []
  },
  {
   "metadata": {},
   "cell_type": "code",
   "execution_count": null,
   "source": "",
   "id": "1087da4a08199b1b",
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T02:54:40.456214Z",
     "start_time": "2025-08-09T02:54:40.435348Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# printre([0.01646966, 0.03255624, 0.06004486, 0.07621205, 0.97623144],[0.85945502, 0.89168254, 0.39898241])\n",
    "# 预测值22.987726671540557,真实值18.193\n",
    "# 预测值25.30597950819623,真实值9.12\n",
    "# 预测值28.03341257118123,真实值9.276\n",
    "# 预测值15.312223386779618,真实值30.243\n",
    "# 预测值16.621474421094778,真实值23.856\n",
    "# 预测值18.469303184319546,真实值35.673\n",
    "# 预测值9.989518771611815,真实值6.586\n",
    "# 预测值11.097051579605028,真实值6.586\n",
    "# 预测值12.257962774797742,真实值8.323\n",
    "printre([0.02299344, 0.16878233, 0.15311101, 0.565999  , 0.36391367],[0.86799529, 0.81583497, 0.05027176])"
   ],
   "id": "13f305ba23439318",
   "execution_count": 77,
   "outputs": []
  },
  {
   "metadata": {},
   "cell_type": "code",
   "execution_count": null,
   "source": "",
   "id": "f275ee343e24d313",
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
